Multi-omic analysis of SDHB-deficient pheochromocytomas and paragangliomas identifies metastasis and treatment-related molecular profiles

Hereditary SDHB-mutant pheochromocytomas (PC) and paragangliomas (PG) are rare tumours with a high propensity to metastasize although their clinical behaviour is unpredictable. To characterize the genomic landscape of these tumours and identify metastasis biomarkers, we performed multi-omic analysis on 94 tumours from 79 patients using seven molecular methods. Sympathetic (chromaffin cell) and parasympathetic (non-chromaffin cell) PCPG had distinct molecular profiles reflecting their cell-of-origin and biochemical profile. TERT and ATRX-alterations were associated with metastatic PCPG and these tumours had an increased mutation load, and distinct transcriptional and telomeric features. Most PCPG had quiet genomes with some rare co-operative driver events observed, including EPAS1/HIF-2α mutations. Two mechanisms of acquired resistance to DNA alkylating chemotherapies were also detected - MGMT overexpression and mismatch repair-deficiency causing hypermutation. Our comprehensive multi-omic analysis of SDHB-mutant PCPG therefore identified features of metastatic disease and treatment response, expanding our understanding of these rare neuroendocrine tumours.


Main
Pheochromocytomas (PC) and paragangliomas (PG) are heritable neuroendocrine neoplasms arising from the sympathetic and parasympathetic nervous system.They are clinically remarkable for catecholamine hypersecretion (e.g.dopamine, noradrenaline, adrenaline) causing morbidity and occasionally death from cardiovascular sequelae 1,2 .PCPG are typically slow growing, but metastases develop in 10-20% of patients 3 .
Current biomarkers are inadequate for predicting metastatic progression, which is currently incurable.
Herein, we conducted an integrated multi-omic analysis of an internationally sourced and clinically wellannotated cohort of SDHB-mutant PCPG including parasympathetic head and neck PG (HN-PG) -a relatively understudied group to date.Importantly, paired primary and metastatic tumours were analysed in a subset of cases.We found novel observations with respect to cell-of-origin, treatment response and clinical outcome and confirmed telomeric dysfunction and other features in metastatic cases.
Review of histopathology and clinical data confirmed disease diagnosis in all cases.Paired synchronous/metachronous primaries were analysed in four cases (E128, E136, E159, E229) and clonal independence confirmed by discordant somatic profiles (Supplementary Figure 1).With respect to clinical behaviour, 40 tumours from 37 patients were ostensibly non-metastatic with at least twelve-months of clinical follow-up (median: 60 months, range 12 -456).An additional eight primary PCPG from eight patients had less than 12 months of clinical follow-up without metastases reported (termed short clinical follow-up).Thirtyfour of 79 patients had confirmed metastatic disease.Twenty-eight primary PCPG were analysed from 26 patients who developed metastases and in six patients at least one paired metastasis was available, confirming the clonal link between the primary and metastatic tumours (termed metastatic primary) (Figure 1C).For 21 primary PCPG from 20 patients, metastasis was reported but metastatic tissue was not available (termed primary -metastasis reported)(Figure 1C).Finally, in eight patients, one or more metastases were analysed but a primary tumour was not analysed.

Molecular profiling confirms genotype-subtype and cell-of-origin transcriptional profile
Molecular subtyping is known to reflect the genotypic features of PCPG 13,[17][18][19][20][21] .To confirm genotypephenotype relationships among A5 tumours we applied UMAP clustering to WTS (n=91), small-RNA-seq and DNA methylation datasets.Published RNA-seq data across a spectrum of PCPG genotypes was included for comparison 13,14,20 .Here, we used the C1/C2 annotation based on seven PCPG gene-expression subtypes previously defined by single-nuclei and bulk gene-expression analysis 22 .As expected, all A5 SDHB-mutant tumours clustered among C1A (SDHx) tumours by WTS.UMAP of the small-RNA-seq data showed a similar clustering pattern; however, four A5 tumours from two patients clustered with an outlier group consisting of unrelated genotypes (Figure 2B).This small-RNA-seq outlier group included C2Bi (MAX) PCPG that frequently exhibit loss-of-heterozygosity of chr14q and silencing of the imprinted DLK1-MEG3 miRNA cluster located on chr14q31-32 20 (Supplementary Figure 3).By DNA methylation profiling, SDHB-mutant PCPG independently of other PCPG genotypes with the exception of one tumour (E229-P1) (Figure 2C), and SDHB-mutant tumours showed genome-wide hypermethylation (Supplementary Figure 4), consistent with a prior report 9 .Unsupervised clustering therefore confirmed the expected molecular profile of the A5 SDHBmutant tumours as well as the pathogenicity of constitutional SDHB variants.
With respect to neurotransmitters, TPH1 (tryptophan hydroxylase), essential for serotonin (5-HT) synthesis, was overexpressed in the HN-PG/non-chromaffin cell tumours (Figure 2G, Supplementary Figure 9).TPH is known to be expressed in glomus cells of the carotid body (presumed the HN-PG cell-of-origin) 29 , and 5-HT storage vesicles have been observed in some HN-PG, although urinary 5-HIAA is not elevated in most HN-PG patients 30 .Meanwhile, while immunohistochemical detection of choline acetyltransferase (CHAT) has been reported in HN-PG 31 , CHAT mRNA was not overexpressed in these tumours (Supplementary Figure 10).

Somatic TERT mutations, structural rearrangement and promoter methylation
Hotspot pTERT mutations (chr5:1295113G>A, chr5:1295135G>A, GRCh38) were identified in 14 metastatic cases (17.7%).Targeted deep sequencing of pTERT excluded potential sub-clonal mutations in the other cases.snATAC-seq data showed a promoter peak central to the de novo ETS binding site 39 in three pTERT mutants tested and all snATAC-seq reads had the somatic variant, confirming open chromatin (Figure 4A).Chromothripsis of chromosome 5 created an AFF4-TERT fusion in one case (E123) (Figure 4B).WGS supported a DNA breakpoint in AFF4 3-prime of untranslated exon 1 and DNA sequence 390 bp upstream of the TERT transcription start site (Figure 4Bii).RNA-seq split-reads supported a fusion transcript involving AFF4 exon 1 spliced to TERT exon 3 (Figure 4Biii).This validates our prior report of TERT structural rearrangements in PCPG 16 but to our knowledge is the first TERT-fusion reported in PCPG.
All TERT-altered cases had elevated TERT expression with the highest expression observed in the AFF4:TERT fusion tumour.TERT copy-number gain did not correlate with TERT overexpression in the absence of TERT mutation or fusion events (Figure 4C).We found TERT promoter hypermethylation including the THOR domain (probes cg11625005, cg10767223) in TERT-altered PCPG, as previously described 15 (Figure 4D).
TERT promoter hypermethylation and concordant overexpression was also observed in a single metastatic TERT wild type tumour (E171-M1) (Figure 4D).
Alternative lengthening of telomeres (ALT) is a known feature of ATRX-altered PCPG and other cancer types 12 .C-circle analysis was done to predict ALT activity 40 and C-circles were detected in all ATRX-altered and one TERT-altered primary tumour (E143-P1) (discussed below)(Figure 5B).As expected, the telomeric tumour:blood ratio, a proxy for tumour telomere length, was higher in ATRX-altered compared to TERT-altered or TERT/ATRX wild type tumours (Student's t-test p-value < 0.00025) with no significant difference in telomeric DNA content observed between the latter two groups.Long telomeres in the absence of C-circles were detected in four ATRX/TERT wild type PCPG, including two HN-PG.
ALT features can include aberrant telomeric variant repeat (TVR) usage, intrachromosomal insertion of telomeric DNA and expression of telomeric repeat-containing RNA (TERRA) 41 .Relative TVR usage was altered in C-circle positive versus negative tumours (Figure 5C).Normalized to total read count only three TVRs (TTAGGG, GTTGGG, GTAGGG) remained significant after false-discovery rate correction (Benjamini & Hochberg, p<0.1).When normalized to the telomeric reads the canonical TVR (TTAGGG) was proportionally more abundant in C-circle positive compared to negative tumours.Only one intrachromosomal telomeric insertion was detected in an ATRX altered case (data not shown), while TERRA expression was elevated in ATRX altered tumours, although was not statistically significant (Two-sided Student's t-test p>0.05)(Figure 5D).
TERT/ATRX alterations are associated with metastatic progression and late somatic events TERT/ATRX alterations were found in 21/34 (61%) metastatic cases and significantly associated with metastasis (Pearson's Chi-squared p-value 9.349e-08).Only one non-metastatic tumour (E156-P1) had a pTERT mutation; however, this case had a short clinical follow-up (Figure 6A).Evidence for late development of pTERT mutations was found in two cases: 1) E146 where the primary was TERT wild type but the paired metastasis was pTERT mutant and 2) case E143 where a subclonal pTERT mutation was detected in the primary (VAF = 10.3%) but higher VAF was detected in two paired metastases (E143-M1 VAF=33%, E143-M2 VAF= 28%), indicating all metastatic cells were pTERT mutant (Figure 6B).Interestingly, the subclonal pTERT mutant primary (E143) also had detectable C-circles and longer telomeres, but these ALT features were absent in both metastases.The data therefore suggests the pTERT mutation and ALT features had coexisted but were likely within different cell subpopulations, while only the pTERT mutant cells metastasized.
Ten primary tumours from patients with metastatic disease did not have TERT/ATRX-alterations, but the late acquisition of TERT/ATRX alterations and presence in metastatic cells could not be excluded given that the metastatic tissue was not sequenced.However, metastasis in the absence of TERT/ATRX alterations was also evident in three metastatic tumours for two patients (E158, E159).These metastases had no other evidence of a TLM (i.e., they did not overexpress TERT, and had neither detectable C-circles nor an increased telomere length ratio) and had no recurrent cancer gene alterations that may explain metastatic progression.Interestingly, we found the tumour mutation burden was significantly elevated in TERT/ATRX-altered PCPG (One-tailed Student's t-test p<0.001) consistent with a previous study 11 , but metastases from the TERT/ATRX wild type cases had a low mutation burden (Figure 6C).Interestingly, while it is also known that a larger primary tumour is associated with metastatic risk, only patients with TERT-altered tumours had larger primaries when contrasting to patients without TERT/ATRX alterations (Two-sided students t-test, p<0.05)(Figure 6D).Finally, survival analysis showed a significant association with poorer outcome in the TERT/ATRX-altered group (Supplementary Table 1, Supplementary Figure 14).

Transcriptional patterns in TERT/ATRX-altered and metastatic PCPG
We next assessed transcriptional changes in TERT and ATRX-altered PCPG, respectively, by contrasting each group independently to non-metastatic PCPG and directly to each other.Furthermore, we contrasted all metastatic and non-metastatic PCPG in a separate analysis (Figure 7A).All HN-PG/non-chromaffin tumours (defined by transcriptional profile) were excluded given the absence of TERT/ATRX alterations in these tumours.Furthermore, we excluded low tumour purity cases (<50% estimated by WGS data) and a subclonal TERT-mutant case.Contrasting the TERT and ATRX-altered tumours to non-metastatic tumours we found 1,152 and 1,448 differentially expressed genes, respectively (Supplementary Data 6).Some of these genes were common to both analyses (n=273), as well as when contrasting all metastatic and non-metastatic tumours (n=180).Consistent with prior studies, genes overexpressed in metastatic PCPG were enriched for cell cycle and proliferation gene ontologies 11,22 .MKI67 was among the top overexpressed genes in keeping with elevated Ki67 IHC staining (Figure 7B).Other genes included cell cycle regulators (e.g.TOP2A, BUB1, FOXM1, AURK2B, CDK1), that were highly correlated with MKI67 expression (r>0.88).Other genes overexpressed in metastatic PCPG included the polycomb repressor complex 2 gene EZH2, as previously described 22 , and the transcription factors OTX1 and TPX1.
Long non-coding RNAs and pseudogenes were enriched in genes differentially expressed in ATRX-altered tumours (Pearson's Chi-squared p-value=9.9e-14)(SupplementaryFigure 15).The serine/threonine protein kinase RIPK4 was also overexpressed, while genes downregulated in ATRX altered tumours included the p53dependent G2M cell cycle regulator RPRM and the GTP-binding gene DRG2 -the latter gene also implicated in G2M cell cycle checkpoint control (Figure 7C) 42 .Other genes down-regulated in ATRX-altered tumours included USE1, SULT4A1 and the nuclear-encoded mitochondrial gene COX17.Many differentially expressed genes were also within differentially DNA methylated regions (Figure 7C).Gene-expression in neoplastic cells was confirmed using the snRNA-seq data and differential expression in ATRX mutant tumours was validated using independent bulk RNA-seq datasets 13,14 (Figure 7D).Meanwhile, few genes were uniquely expressed in TERT-altered tumours but included IRX3, SDK1 and TRIP13 (Figure 7E).

Treatment-related mutagenesis and resistance to DNA alkylating chemotherapy
DNA alkylating chemotherapies (e.g.dacarbazine and temozolomide) are common treatments for metastatic PCPG 43 .Two patients had paired tumours taken before and after cyclophosphamide, vincristine and dacarbazine (CVD) treatment.Case E169 had a para-aortic AT-PGL with distant lung, liver, and bone metastases (Figure 8A).Tumours resected before (lung, E169-M1) and after (spinal disease, E169-M2) 23 cycles of CVD were used for genomic analysis.Paired E169 tumours had a high number of shared SNVs (n=3158) including mutations in TERT and BCOR; however, the number of SNVs doubled in the post-CVD (n=12702) tumour including mutations in TP53, RPL5 and POLE (DNA Polymerase Epsilon) (Figure 8B).SBS signatures also changed in the post-treatment sample with increased representation of SBS8, SBS12, SBS16 and SBS21 as well as indel signature ID1 (Figure 8C).However, elevated SBS11 associated with DNA alkylating chemotherapy and SBS14 associated with DNA polymerase epsilon-deficiency was absent.
Importantly, MGMT was overexpressed in the post-treatment tumour (E169-M2) (Figure 8D).MGMT overexpression is a known temozolomide resistance mechanism in glioblastoma cells 44 and therefore potentially an acquired dacarbazine resistance mechanism in PCPG.
A second case (E167) had a left adrenal primary PC and then 14 years later developed HN-PGL and metastatic disease in the hip (Figure 8E).The patient progressed with development of liver metastases after 131 I-MIBG therapy and was subsequently treated with 22 cycles of CVD.Further progression occurred with new intracranial and osseous disease.Liver and intracranial metastases represented disease pre-and post-CVD treatment respectively, with clonal relatedness confirmed by common SNVs (n=2407) including a pTERT mutation (Figure 8F).MGMT expression was low in both tumours (Figure 8D).The post-CVD metastasis acquired more than a million additional SNVs corresponding to a strong SBS11 signature.Interestingly, the post-treatment intracranial metastasis (E167-M2) had a somatic NRAS (NP_002515.1:p.G12A) mutation and predicted to have elevated MAP kinase pathway activity by WTS (Supplementary Figure 16).Importantly, a protein truncating MLH1 splice-site mutation (chr3:g.37091976G>A) was detected in the post-CVD metastasis with concordant loss-of-heterozygosity and low MLH1 expression (Figure 8F).Canonical SBS mismatch-repair (MMR) mutational signatures (e.g.SBS6, SBS14, SBS21) were not detected in the post-CVD treated tumour; however, insertion/deletion signature ID2 was elevated and insertion/deletion events at microsatellite sites revealed a high MSI-score (11.7 microsatellite indels per mega base) in the post-CVD tumour (E167-M2), which is above the prescribed microsatellite instability threshold (>4)(Supplementary Figure 17).Pathogenic variants in MMR genes (MLH1, MSH2, MSH6, PMS2) as well as treatment-related hypermutation and SBS11 have been described in temozolomide-treated glioblastomas 45,46 and more recently also pancreatic neuroendocrine tumours 47 .To search for evidence of MMR-deficiency in additional PCPG we analysed variant data from AACR Project GENIE (release 13), representing comprehensive cancer panel data for 75 PCPG.We identified one hypermutated case with a dominant C>T mutational pattern and an MSH6 variant predicted to be a driver by Cancer Genome Interpreter (GENIE-DFCI-001077-11300, NP_000170.1:p.G162E) (Figure 8I-J).Furthermore, we performed 523-gene panel sequencing of plasma cellfree (cf) DNA from a patient with SDHB-associated PCPG with high disease volume as shown by 18 F-FDG-PET (Figure 8K).Blood was taken after three cycles of temozolomide combined with 177 Lutetium-DOTAoctreotate.The patient had high fraction of circulating tumour (ct) DNA (45%) (Figure 8L).Hypermutation was evident in cfDNA with a predominance of C>T transitions corresponding to SBS11 (Figure 8J), although most C>T mutations had low VAF (~0.5%) in contrast to truncal mutations involving pTERT, MYCN, NUTM1 and MAP2K1 (VAF>10%) (Figure 8M).A mutation in an MMR gene was not detected; however, the hypermutated clone likely represented a small fraction of total disease volume therefore detection of an MMR gene mutation may have fallen below the detection limit of the cfDNA assay (VAF>0.5%).

Discussion
We performed a comprehensive genomic analysis of SDHB-associated PCPG thereby creating a data resource to understand development of metastatic disease in these patients.Consistent with prior studies, most PCPG had stable genomes with few recurrently mutated cancer genes 11,13,14,20 .Importantly, we confirm SDHBmutant PCPG are molecularly distinguishable from other PCPG genotypes, while parasympathetic (nonchromaffin) tumours including HN-PG had a very distinct molecular profile.Our data confirms TERT and ATRX alterations are the most common secondary driver events in PCPG and are enriched in metastatic disease 11,13,15 .However, a small number of metastatic PCPG lack TERT/ATRX alterations, suggesting yet unidentified mechanisms of metastatic progression in these cases.Mutual exclusivity of TERT/ATRX-alterations in PCPG suggests these events have a common and potentially redundant role in metastasis.Both genes are involved in telomere maintenance -a hallmark of cancer 48 .
However, given the majority of PCPG, including a subset of metastatic PCPG, have no obvious telomere lengthening mechanism (TLM), how, or indeed whether, these tumours maintain their telomeres is unknown.
While methods such as the Telomere Repeat Amplification Protocol (TRAP) 49 and the C-circle assay are required to exclude the presence of a known TLM in tumour cells, absence of TLM has also been reported in other neuroepithelial cancers 50,51 .Meanwhile, while ATRX/TERT alterations were strongly associated with metastatic progression, the presence of long telomeres in the absence of TLM appears unlikely to be sufficient to drive metastatic progression 15 .
ATRX is a SWI-SNF chromatin remodelling protein with multiple functions beyond telomeric regulation, while TERT is also thought to have non-canonical functions 52 .We observed that increased cell cycle activity appears to be common to ATRX/TERT altered PCPG, which is consistent with prior studies 11,22 and in keeping with increased Ki67 staining within metastatic PCPG 53 .Experimental models implicate ATRX and TERT in regulating cell cycle activity.For instance, ATRX binds to regulatory elements of cell cycle genes including CHEK1, as previously reported in a glioblastoma model 54 .Intriguingly, we found repression of RPRM and DRG2 in ATRX altered PCPG, which are known to regulate the G2M cell cycle checkpoint 42,55 .It remains to be determined if RPRM and DRG2 are ATRX target genes.TERT can also promote cellular proliferation and has a pro-survival function.For instance, cellular proliferation can be inhibited by wild-type reversion of a pTERT mutation in a glioblastoma cell line 56 .Furthermore, catalytically inactive TERT can rescue TERT-null cells, indicating a pro-survival function that is independent of telomerase activity 57 .How TERT promotes cellular proliferation is poorly understood and may be cell type and context dependent.Prior studies have shown TERT occupancy at gene-promoters including tRNAs 58 while TERT has also been shown to interact with oncogenic pathways involving MYC, WNT/B-catenin or NFkB 52 .Given few differentially expressed genes were restricted to TERT-altered PCPG, experimental models will be required to understand any noncanonical TERT functions in PCPG.
ATRX/TERT alterations have potential implications for treatment and surveillance.An increased sensitivity of ATRX-altered tumour cells to DNA damage response (DDR) inhibitors targeting ATM 54 , WEE1 59 , ATR and PARP 60,61 has been reported in other neuroectodermal cancers, and therefore might be effective in PCPG.
Sensitivity to radiosensitizing DDR-targeting drugs in combination with radionuclide therapies (e.g. 177Lutetium-DOTA-octreotate or 131 I-MIBG) also warrants investigation in the context of ATRX/TERT alterations 62,63 .Although we confirmed TERT/ATRX-alterations are associated with metastatic disease, the utility of these mutations as predictive biomarkers may be limited by the late timing of these events.Detection of TERT/ATRX alterations by cfDNA analysis may be a viable alternative to tissue biopsies, although the broader utility of cfDNA detection in PCPG still requires further validation.
DNA alkylating agents are commonly used for PCPG but there are currently no validated biomarkers to predict treatment response.We identified hypermutation together with MMR mutations (e.g.MLH1, MSH6) or MGMT overexpression in dacarbazine or temozolomide-treated PCPG.These observations are potentially clinically important.Firstly, monitoring PCPG patients for development of MMR mutations, SBS11 or potentially MGMT overexpression may indicate development of treatment-resistant cells.Secondly, MMRmutations may sensitize tumour cells to other treatments.For instance, hypermutation can theoretically generate neoantigens causing T-cell recognition in tumours, making them more immunogenic and potentially more sensitive to immune checkpoint inhibition (ICI), although this theory is tempered by observations in glioblastoma where temozolomide-induced hypermutation did not translate to ICI-response 46 .Synthetic lethality could be another approach in MMR-deficient tumours, including use of WRN helicase inhibitors 64 .
Further investigation is required to determine whether these treatments are efficacious in MMR-deficient PCPG.
Although therapeutically actionable mutations are rare in SDHB-mutant PCPG, co-operative drivers can still represent potential treatment opportunities in some patients.In particular, somatic EPAS1 variants in SDHBmutant tumours may have translational relevance.It is currently unknown whether HIF2α inhibitors (e.g.belzutifan) are effective for SDHB-mutant PCPG, but it is plausible that co-operative EPAS1 mutations could make these tumours more sensitive 65 .Screening for co-operative somatic EPAS1 mutations in SDHB-mutation carriers should therefore be employed for future HIF2α inhibitor clinical trials.Although we have reported the most comprehensive analysis of SDHB-mutant PCPG to date, future studies will be needed to validate our findings and identify additional low frequency genomic features that may be used as clinical biomarkers and future treatment targets.

Tissue sample review and nucleic acid extraction
Tissue samples were assessed for tumour cell content by a single pathologist (AG) by haematoxylin and eosin staining of 4uM cryosections.Ten serial 10uM sections were taken for nucleic acid extraction adjacent to the H&E reviewed sections.If required, microdissection of tissue was performed on ethanol fixed and stained fresh tissue sections to dissect tumour regions from the sectioned material.DNA was extracted from tissue samples and matched patient whole blood using the DNeasy® Blood & Tissue kit (Qiagen, Germany) according to manufacturer's protocol.Total RNA was extracted from tissue using the miRNeasy mini kit (Qiagen, Germany) according to the manufacturer's protocol.Nucleic acids were quality control checked by Agilent TapeStation microfluidic electrophoresis (Agilent, USA) and NanoDrop spectrophotometer (Thermo Fisher Scientific, USA).

Targeted Amplicon Sequencing
To validate mutations of interest, targeted amplicon sequencing was done using custom primers compatible with universal indexed primers for Illumina cluster generation.PCR amplification for the region of interest was done in a single well (for each patient sample) and the expected PCR amplicon product then verified by agarose gel electrophoresis.Multiple amplicons for each patient sample were pooled, purified, and then barcoded with the index primers, quantified by fluorometer (Qubit dsDNA High-sensitivity, USA) and then sequenced on a MiSeq using a V3 300 cycle kit (Illumina, USA).FASTQ files were aligned to human genome (GRCh38) using BWA-MEM and BAM files were visualized in IGV 66 to confirm variant of interest and record the variant allele frequency.

Library preparation
Libraries were prepared using the Illumina® TruSeq™ DNA Nano library preparation method according to the manufacturer's instructions at The University of Melbourne Centre for Cancer Research (UMCCR) using 200ng input DNA and a 550 base pair insert size.Samples were sequenced in separate batches on the Illumina® Nova-Seq 6000 according to manufacturer's instructions (Illumina, USA).

WGS alignment and small variant calling
BCL files were demultiplexed and converted to FASTQ files with Illumina® bcl2fastq (version 2.20.0.422) with default settings (including adapter trimming).Alignment and variant calling steps were run using the bcbio-nextgen cancer somatic variant calling pipeline (version 1.2.4-76d5c4ba) (https://github.com/bcbio/bcbio-nextgen). The versions of all programs used by bcbio-nextgen are additionally shown in Supplementary Table 2. Briefly, the tumour and normal sequencing data in FASTQ format were first processed by Atropos (version 1.1.28) 67to clip homopolymers (minimum 8 bases) from the ends of reads.Trimmed FASTQ files were aligned to the human genome (GRCh38) with BWA-mem (version 0.7.17) 68 with predominantly default settings.The two exceptions were that the -M flag was enabled to mark shorter split read hits as a secondary alignment and seeds with > 250 occurrences were skipped (reduced from the default 500).
Germline variants were called by vardict (version 1.8.2) 69 , GATK HaplotypeCaller (GATK version 4.1.8.1) 70,71 and Strelka2 (version 2.9.10) 72 callers.Somatic variants were called by vardict, Mutect2 (see GATK version) 73 and Strelka2 (version 2.9.10) variant callers.In both instances a variant was accepted as true if it was detected by 2 of the 3 callers.Somatic variants were further post-processed by umccrise (version 1.2.0-rc.1,https://github.com/umccr/umccrise).To further reduce spurious variant calls, a variant blacklist was generated by counting supporting reads for all detected somatic variants across all normal BAM files using SAMtools 74 (version 1.12) mpileup.Any somatic variants that were supported by three or more reads in more than two normal controls were excluded.

WGS copy number and structural variant calling
Structural variants (SV) were called with GRIDSS (version 2.8.0) 75 .GRIDSS output was annotated for repeat sequences using RepeatMasker (version 4.1.5)and filtered for somatic SVs using GRIPPS (version 1.11).
To ensure high quality somatic calls, all SVs passing GRIPSS filtering were examined in the raw output of GRIDSS across all samples.Any SV (based on exact genomic coordinates) that was observed (whether ultimately rejected or passed) in more than two unrelated samples was removed.Finally, SVs were required to have either a minimum of three supporting split reads or paired reads, or a combination of at least one split and one paired read.Somatic copy number profiles were generated using PURPLE (version 3.1, https://github.com/hartwigmedical/hmftools)from input generated by COBALT (version 1.11), AMBER (version 3.5), and GRIDSS.Structural variants were further annotated with Linx (version 1.16) using the output from PURPLE.

Telomere length estimation and telomere variant repeat detection
Telomere content was estimated from WGS for each tumour and normal pair using TelomereHunter (version 1.1.0) 76with default parameters.Telomere content was computed from the number of intratelomeric reads normalized by the total number of reads with a GC composition like that of telomeres.Telomeric enrichment in tumours was defined as a log2 ratio of tumour telomere content to blood telomere content greater than 0.5.TelomereHunter was also applied to WTS for each tumour to estimate TERRA content.

Mutational signature analysis
For each sample, the single base substitution (SBS), doublet base substitution (DBS), and small insertions and deletions (ID) mutation count matrices were extracted from variant call format (VCF) files containing somatic variant calls using the MutationalPatterns (version 3.10.0) 77package for R. The optimal linear combination of the COSMIC V3 78 SBS, DBS, and ID mutational signatures that best reconstructed their respective profiles was then estimated using the fit_to_signatures function.

Microsatellite instability assessment
Insertions and deletions in microsatellite regions were assessed using PURPLE (version 3.1).Samples with more than four microsatellite Indels per megabase were considered to have microsatellite instability as per the software default setting.

Identification of recurrent copy number alterations
The Genomic Identification of Significant Targets in Cancer (GISTIC, version 2.0) 79 algorithm was used to identify recurrent copy alterations.The input data for GISTIC was modified from the somatic segmentation output from PURPLE by taking log2 of the segmental copy number incremented by a small offset (0.01) and subtracting one to center the diploid state around zero (log2(CN + 0.01)-1).GISTIC was run in broad mode with segments containing less than 50 markers undergoing joining.The broad length cutoff was set to 75% of a chromosome and arm-peel was activated.The analysis was applied independently to sample subsets including HN-PG/parasympathetic PG, non-metastatic and metastatic sympathetic PCPG.In cases where a patient had more than one primary tumour, one was chosen at random, in cases where a paired primary and metastasis were present, the metastasis was used.

Identification of driver genes
The dNdScv R package 33 was used to assess somatic mutations to detect genes under positive selection.
Before running the analysis, quality filtered somatic variant calls were lifted over from the GRCh38 to HG19 genome build using the rtracklayer package for R with a liftover chain obtained from UCSC (hg38ToHg19.over.chain).The resulting mutation set was analysed with dndscv function using default parameters.Additionally, the HG19 mutation set was used to run the OncodriveFML 34 driver detection tool via the web-based application interface (https://bbglab.irbbarcelona.org/oncodrivefml/home,accessed 5 th February 2024).The algorithm was run using genomic region subsets encompassing coding sequences or promotor regions with Combined Annotation Dependent Depletion (CADD) scores version 1.0.

Library preparation
Samples were prepared using the NEB-Next directional RNA-Seq kit (NEB, USA) and underwent 150bp paired-end sequencing on the Illumina NovaSeq 6000 (Illumina, USA) according to manufacturer's instructions.

Read quantification
RNA sequencing reads were aligned with STAR 80 (version 2.7.9a) using a two pass approach.First, all data was aligned using a STAR reference built with the Gencode human reference genome with decoy HLA sequences (GRCh38) and Gencode gene transcript annotation file (version 36, primary assembly annotation).Following alignment, the splice junction output from all samples was merged.Supporting read evidence for each junction was summed and junctions with less than 10 unique supporting reads were removed.The resulting junction file was used to build a new STAR reference and the alignment was repeated.Gene level feature counts were obtained with HTSeq-count (version 0.11.3) using the previously described Gencode transcript annotation file.All parameters were default except mode which was set to intersection-nonempty.

Differential gene expression
To contrast gene expression between parasympathetic and sympathetic (including pheochromocytoma) paraganglioma, samples were assigned to one of two groups based on clinical annotation and UMAP clustering.Samples where UMAP clustering did not reflect the clinical annotation or tumour purity was below 50% were excluded from the contrast.To perform differential expression analysis, transcript counts obtained from HTSeq-count were converted to a digital gene expression object using the edgeR (version 3.42.4) 81package for R. Lowly expressed genes were removed using the filterByExpr function with default parameters.Library scaling factors were then computed using the calcNormFactors function.A design matrix was constructed without an intercept to include the differential group and patient gender as a covariate.Normalization factors were applied and correction weights for count-dependent variance were computed using the voom function from the limma 82 (version 3.56.2) package for R. To account for multiple samples from the same patient, the inter-duplicate correlation was computed using the duplicateCorrelation function with the patient identifier as the blocking factor.Linear model fitting was performed using the lmFit function supplied with the duplicate correlation computed above and the patient identifier as the blocking variable.Contrasts were fitted using the contrasts.fitfunction and test-statistics were computed using the eBayes function.
To contrast metastatic, TERT-altered, and ATRX-altered tumours with non-metastatic primary tumours samples were assigned to groups based on clinical behaviour.The clinical behaviour groups were non-metastatic primary tumour, primary tumour where metastatic disease was present, metastatic lesion, or a tumour with short clinical follow up.These groups were then further segregated based on the presence or absence of a TERT or ATRX mutation.Parasympathetic tumours and tumours with low purity (<50%) were excluded from all the contrasts.Samples with an ATRX/TERT alterations where the allele fraction was less than 25% were excluded from genotype-centric contrasts, and samples with short clinical follow-up were excluded from contrasts involving clinical behaviour.See Supplementary Data 2 for a full list of contrast participation.Differential expression was performed as described above.

RNA fusion detection
Whole transcriptome sequencing data aligned with STAR aligner (described above) was analysed with the Arriba (version 2.1.0)gene fusion detection tool using the supplied blacklist, known fusion, and protein domain files with default parameters.

Harmonization with published RNA-seq datasets
Read counts previously published by Fishbein et al. were obtained from the Genomic Data Commons using the TCGAbiolinks 83 (version 2.29.1) package for R with these parameters: project = "TCGA-PCPG", data.category= "Transcriptome Profiling", data.type= "Gene Expression Quantification", workflow.type= "STAR -Counts".Additional data previously published by Flynn et al 14 was obtained from the European Genome-Phenome Archive (EGAS00001005861).Raw counts from the publicly available data were merged with our dataset retaining genes present in all three datasets.Count matrices were converted to a digital gene expression object with edgeR and lowly expressed genes were removed using the filterByExpr function with default parameters.Library scaling factors were then computed using the calcNormFactors function and counts were transformed to log counts per million using the cpm function.Possible technical variance was ameliorated using the removeBatchEffect function from the limma package where the dataset identifier was supplied as the batch and the primary driver genotype (e.g RET, SHDx, VHL) as the factor to preserve.

Library preparation
Samples were prepared using the NEXTFLEX® Small RNA-Seq Kit v3 (Bioo Scientific).Samples were sequenced at the molecular genomics core (Peter MacCallum Cancer Centre) using 50bp single end sequencing on the Illumina NextSeq 500 (Illumina, USA).

Read quantification
FASTQ files were pre-processed with two rounds of cutadapt (version 1.6.3) to first remove adapter sequences (-a TGGAATTCTCGGGTGCCAAGG, --minimum-length 23) and then remove the first and last four bases (-u 4, -u -4).Following trimming, sequences were aligned to GRCh38 with subread (version 1.6.3)and quantified using featurecounts (version 1.6.3) in stranded mode against the miRBase (version 20) miRNA transcripts 84 .Following read quantification, read-count processing was performed in an identical manner to that described under RNA-Seq analysis.

Harmonization with published small-RNA-seq datasets
Small-RNA sequencing data previously published by Castro-Vega et al were obtained from ArrayExpress (E-MTAB-2833, https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2833).Additional data published by Fishbein et al. were downloaded from the Genomic Data Commons using the GDC Data Transfer Tool in BAM format and converted to FASTQ format using SAMtools (1.10).Alignment and read quantification was performed as described above under 'Read quantification'.Read count harmonization was performed as described for RNA-seq dataset under 'Harmonization with published RNA-seq datasets'.

Methylation array data pre-processing and normalization
The minfi (version 1.46.0)package for R was used to read in a targets manifest and red/green channel IDAT files using the read.metharray.sheet and read.metharray.expfunctions, respectively.Probe detection pvalues were determined with the detectionP function.Raw green/red values were converted to methylation signal and functionally normalized using the preprocessFunnorm function from minfi.Any probes with a detection p-value below 0.01 in any sample or those that fell on a SNP position (using dropLociWithSnps) were removed from the dataset.Furthermore, the maxprobes package for R was used to exclude any crossreactive probes.Illumina probe annotation version 1.0 with GRCh38 coordinate information was obtained from: https://github.com/achilleasNP/IlluminaHumanMethylationEPICanno.ilm10b5.hg38

Probe level differential expression
Probe level differential methylation analysis was performed using the Remove Unwanted Variation workflow with the RUV package via the MissMethyl package for R 85,86 .Firstly, functionally normalized methylation values were converted to M-values using the getM function from the minfi package and negative control probes were obtained using the MissMethyl's getINCs function.For each contrast, two rounds of variation adjustment were performed.In the first round, an RUV model was computed with RUVinverse algorithm using the RUVfit function applied using the Illumina negative control probes as controls.
The rescaled variances were then applied using the RUVadj function and p-values for differential methylation of each probe were obtained with the topRUV function.Probes with a p-value greater than 0.5 were selected for use as controls in the second round of variation adjustment using the RUVfit function and the resulting model was applied using the RUVadj function.

Differentially methylated regions
Identification of differentially methylated regions was performed with the DMRcate 87 package (version 2.14.1) for R. Firstly, functionally normalized methylation values were converted to M-values using the getM function from the minfi package.Contrast and design matrices were generated using the same process described for differential gene expression.For each contrast, individual probes were then annotated using the cpg.annotate function with analysis type set to differential and the default annotation set for Illumina EPIC arrays (ilm10b2.hg19).Differentially methylated regions were then identified using the dmrcate function with a lambda of 1000 and bandwidth scaling factor of 2. Significant DMRs were those with a Fisher multiple comparison statistic below 0.05.

Harmonization with previously published methylation array datasets
Illumina Infinium HumanMethylation450 methylation array data previously made available by Fishbein et al. were obtained from the Genomic Data Commons using the TCGAbiolinks 83 (version 2.29.1) package for R with these parameters: project = "TCGA-PCPG", data.category= "DNA Methylation", data.format= "IDAT", legacy = False.The array data were processed as described under 'Methylation array data preprocessing and normalization' with the exception of the annotation file, which was substituted with the supplied Illumina Infinium HumanMethylation450 annotation.Additional data made available by Letouzé et al. 9 were obtained from the Gene Expression Omnibus using the GEOquery package for R (accession GSE39198 and GSE43293).The GSE39198/GSE43293 data was only available as beta values and was not reprocessed.The beta-values from all three datasets were converted to m-values and merged, retaining only probes available on the Illumina Infinium HumanMethylation27 array design.A design matrix was generated with the model.matrixfunction from the limma package using the primary driver genotype (e.g.RET, SHDx, VHL, Unknown).Batch effects were then adjusted using the removeBatchEffect function from limma where the dataset was supplied as the first batch annotation and patient gender as the second batch annotation.

Data processing and analysis
Additional snRNA-seq data for two normal adrenal medulla samples (E240, E243) previously published by Zethovan et al. 22 were obtained from the European-genome-phenome archive (accession EGAS00001005861).Single-nuclei RNA sequencing reads in FASTQ format were aligned and counted using CellRanger (version 7.2.0).Read-one data for libraries processed with the 3' protocol were restricted to the first 26 bases using the r1-length flag.Feature counting was performed using the 10x supplied GRCh38 transcriptome (version refdata-gex-GRCh38-2020-A) with reads aligning to introns contributing to feature counts.Single-Cell Remover of Doublets (SCRUBLET, version 0.2.3) was used to identify gelbeads containing more than one cell.CellRanger count matrices were loaded using the Read10X and CreateSeuratObject functions from the Seurat (version 4.3.0.1) package for R in conjunction with the SCRUBLET output.For each sample, features present in less than three cells and cells containing less than 500 features or 750 unique RNA molecules, or more than 10 percent mitochondrial RNA were excluded.

Single nuclei ATAC-sequencing
Library preparation with 10x Chromium snATAC-seq was conducted using the "Van Helsing" protocol (dx.doi.org/10.17504/protocols.io.bw52pg8e).Initially, nuclei were isolated following the steps outlined in the previously mentioned "Frankenstein" protocol.Sorted nuclei were collected into a round-bottom 96-well plate containing 100 µL of ice-cold Wash Buffer.The sorted nuclei were then transferred to 0.2 mL PCR tubes; 50 µL of ATAC Wash Buffer-Dig was added to the well to ensure transfer of any remaining nuclei.
The nuclei were then centrifuged at 500xg for 5 minutes at 4°C.The supernatant was carefully decanted, leaving behind approximately 10 µL.Without resuspending the pellet, 100 µL of ice-cold Diluted Nuclei Buffer was added.The sample was centrifuged again under the same conditions, and 100 µL of the supernatant was gently removed in two steps-first 90 µL, then the remaining 10 µL.Another 100 µL of Diluted Nuclei Buffer was added gently.This washing step was repeated, removing the supernatant in two steps to leave behind approximately 7-10 µL.The nuclei were resuspended in this small volume, ensuring thorough mixing by carefully washing the walls of the tube.For quantification, 1-2 µL of the nuclei suspension was diluted 1:5 with Diluted Nuclei Buffer and mixed 1:1 with Trypan Blue.Nuclei were counted using a Countess II FL Automated Cell Counter or a hematocytometer to estimate the expected number of nuclei based on the recovery factor.Finally, 5 µL of the nuclei suspension in Diluted Nuclei Buffer were used directly in the Chromium Single Cell ATAC Reagent Kits protocol (CG000168 Rev A), following the specific recommendations for low input samples to optimize the transposition reaction and subsequent library preparation.. Once processed, snATAC-seq libraries were sequenced on the Illumina NextSeq 500 (Illumina, USA) using 50bp paired-end, duel-indexed sequencing.

Data preprocessing
Raw sequencing data was demultiplexed using cellranger-atac mkfastq (V2.0.0).Cellranger-atac count was used to combine replicate sequencing runs, align single nuclei ATAC reads to the GRCh38 reference genome (GRCh38 2020 A 2.0.0, 10x Genomics), and call initial peak sets for individual samples.
Cellranger-atac aggr was then used to combine samples into a single peak-barcode matrix.

Quality control
Signac 88 (version 1.12.9004) and Seurat (version 4.3.0.1) were used for quality control and downstream analysis of snATAC-seq data.Filters were applied to QC metrics that were independent of peaks, to avoid bias toward under-represented cell types.Transcription start site enrichment scores were calculated using the Signac TSSEnrichment function.Outlier nuclei with over 5000 unique reads were assumed to be doublets or nuclei clumps and removed.Nuclei with TSS enrichment scores below 2 were removed due to having a low signal-to-noise ratio.

Cell type identification
Cell type identification was performed via label transfer in Signac.A gene-by-cell accessibility matrix was computed by calculating the sum of reads aligning to the gene body and promoter regions of protein-coding genes for each cell.Log normalisation was performed on the gene activities with Seurat NormalizeData.
Integration was performed to link snATAC-seq nuclei clusters to defined cell types within the snRNA-seq dataset generated from matched and similar samples 89 .Integration was performed with the FindTransferAnchors Seurat function.Briefly, Canonical Correlation Analysis (CCA) was used to identify shared sources of variation between the snATAC gene activity and snRNA gene expression, and project the two datasets into the same reduced-dimensional space.Mutual nearest neighbours (cell pairs from each dataset within each other's neighbourhood) were identified and used as integration anchors between the two datasets.These anchors were then used to transfer cell-type labels from the snRNA reference dataset to the to snATAC-seq dataset, using the Seurat TransferData function, enabling prediction of snATAC-seq nuclei cell types.Cell type identification was confirmed by gene activity of canonical cell-type-specific marker genes.

Peak calling and visualisation
Open genomic regions were defined using pooled reads for each cell type and tumour sample with > 100 total nuclei individually using the MACS2 90 peak caller (version 2.2.9.1) implemented in the Signac CallPeaks function.Resultant fragment pileup bedgraph files for each cell type were visualized using the SparK command line tool 91 .

Cell-Free DNA panel sequencing Sample collection
Patient blood was collected in Cell Free DNA BCT CE tubes (catalogue number 218996, Streck, USA).
Within 2 hours of collection, whole blood was centrifuged at 1900g (brake off) for ten minutes to separate the plasma.Plasma was aspirated and aliquots were centrifuged at 20000g to pellet residual cells or debris.
The plasma was stored at -80°C until DNA extraction.Cell free DNA (cfDNA) was extracted from 5mL of plasma using the Qiagen QIAamp circulating nucleic acid kit (catalogue number 55114, Qiagen, The Netherlands) according to manufacturer's instructions.

Library preparation
A sequencing library was generated with 30ng of input cfDNA using the Illumina TruSightOncology 500 ctDNA library preparation kit (catalogue number 20039252, Illumina, USA).Input quantity was determined based on the mononucleosome fraction (75-300bp) using the Agilent TapeStation4200.Library preparation was performed as per manufacturer's instructions with the addition of a 1.2x bead clean-up of the amplified enriched library.The fragment distribution of the enriched library was assessed on the TapeStation4200, and library quantification was performed using the KAPA Illumina library quantification kit (catalogue number 07960140001, Roche, Switzerland) and analysed with QuantStudio6 (ThermoFisher, USA).

Sequencing and data processing
Sequencing was performed at the University of Melbourne Centre for Cancer Research Centre using the Illumina NovaSeq 6000 (Illumina, USA) to achieve a depth of coverage greater than 35000x.Data was analysed with Illumina DRAGEN TruSight Oncology 500 ctDNA Analysis Software v1.2 on the Illumina Connected Analytics platform (ICA v1).PierianDx Clinical Genomics Workspace (v6.27.0.1) was used for variant interpretation.

C-Circle Assay
Isothermal rolling circle amplification and detection of C-circles was performed as described 40,94 .For each tumour sample, 20 ng DNA was incubated with 10 μl 0.2 mg/ml bovine serum albumin, 0.1% Tween, 4 mM DTT, 1 mM each dATP, dGTP, dCTP and dTTP, 1× Φ29 Buffer , with or without Φ29 polymerase 7.5 U (New England Biolabs) at 30 °C for 8 hours, then at 65 °C for 20 min.Reaction products were diluted to 100 μl with 2× SSC and dot-blotted onto a 2× SSC-soaked Biodyne B nylon membrane (Pall Corporation).The membrane was air-dried, twice crosslinked using 1,200 J UV (Stratagene) and hybridized overnight at 37°C with end-labelled 32P-(CCCTAA)3 and PerfectHyb Plus hybridization buffer (Sigma).The membrane was then washed at 37°C three times in 0.5x SSC/0.1% SDS and exposed to a storage phosphor screen, which was imaged using an AmershamTM TyphoonTM laser scanner (Cytiva) and analysed by densitometry using ImageQuant software (Molecular Dynamics).C-Circle levels for each tumour sample were reported as the difference in the dot-blot intensities resulting from reaction of the sample's DNA with and without Φ29 polymerase.

Statistical analysis
Student's t-tests and Pearson's Chi-squared tests for count data were performed using the t.test and chisq.testfunctions, respectively, from the base stats package for R (version 4.3.1)telomere content ratio (log2) greater than (red) or less than (grey) 0.5 (D) Reads containing telomeric sequences were counted from WTS data using TelomereHunter.The telomeric content (y-axis) was computed as the number of unmapped reads containing telomeric sequence times 1,000,000 divided by the total number of reads with a GC content similar to telomeric repeats.Data point colours, x-axis, and statistical testing are as described in panel B. The lower and upper hinges of each boxplot correspond to the first and third quartiles, respectively, and the median value is marked.The whiskers extend to the largest and smallest value no greater than 1.5 times the interquartile range above or below the upper and lower hinges, respectively.The dimensions (centimetres, y-axis) of the largest primary tumour reported for each patient.Patients are stratified by the presence of a TERT/ATRX mutation and the presence (M) or absence (NM) of metastatic disease.The lower and upper hinges of each boxplot correspond to the first and third quartiles, respectively, and the median value is marked.The whiskers extend to the largest and smallest value no greater than 1.5 times the interquartile range above or below the upper and lower hinges, respectively.A one-tailed Student's t-test used to test for significance.
Tumour DNA samples were normalized to approximately 200ng of DNA in 45µL.Bisulfite conversion was performed with the EZ-96 DNA Methylation MagPrep (Zymo Research, California, USA).Bisulfite converted DNA was processed on the Illumina Human Methylation EPIC BeadChip (Illumina, California, USA) following the Illumina Infinium HD Methylation assay protocol (AGRF NATA scope GGTMN00263)

Figure 1 :
Figure 1: Multi-omic analysis of SDHB-mutated PCPG (A) Outer ring: Overview of the analytical methods applied to the cohort annotated with the number of tumours analysed.Inner ring: The number of tumours and patients (in parenthesis) analysed from each anatomical location.Where a metastasis was analysed, the location of the primary tumour is indicated.(B) A summary of the combination of analytical methods applied to each tumour with respect to anatomical primary location.The upper panel indicates the total number of samples from each anatomical location analysed with the respective combination of assays (lower panel).(C) Anatomical location (dot colour) and clinical behaviour (tile colour) of each tumour analysed.Paired samples from the same patient are joined by a line.(D) Summary of the germline SDHB mutations detected by WGS across the cohort and their respective position within the protein amino acid sequence.(i) Amino acid changes for single nucleotide variants and small insertion/deletion events (NP_002991.2). (ii) The total number of patients observed with each amino acid change where bar colour indicates the clinical disease course.(iii) Schematic of the consequence of splice donor mutations.(iv) Regions of the SDHB protein sequence deleted by large scale structural events.Bar colour indicates the clinical disease course and the number of patients affected is indicated within the bar.

Figure 2 :
Figure 2: Genomic profiling separates sympathetic and parasympathetic PGL UMAP dimensional reduction was used to cluster WTS (A), small-RNA-seq (B) and DNA methylation data (C) with previously published data.UMAP clustering was also repeated for respective A5 data types in isolation (D-F).(G) Differential expression profiling was performed between abdominal-thoracic PCPG and HN-PG.The heatmap shows CPM (log2, z-score) expression values for each tumour (x-axis) for the top differentially expressed genes (adjusted p-value < 0.05, ranked by log-fold-change, top and bottom 30 genes are shown)(y-axis).The top annotation bars indicate the suspected cell-of-origin based on UMAP clustering and the anatomical location of the tumour, respectively.The left annotation bar indicates whether the gene was also in a differentially methylated region for the same contrast.(H/I) Spatial distribution of adjusted p-values (y-axis, -log10) along chromosomes 7, 12, and 17 (x-axis) from differential expression (H) and probe-level differential methylation analysis (I) between sympathetic PCPG and HN-PG.(J) Expression of catecholamine biosynthesis and processing pathway genes.Line colour indicates the anatomical location of the tumour and sub-panels segregate tumours based on which catecholamines were above upper normal limit during clinical testing.Expression data for the A5 cohort were combined with a larger compendium of publicly available PCPG representing the different PCPG subtypes before values were normalized to Z-scores.Only A5 tumours are shown.

Figure 3 :
Figure 3: Genomic and clinical features of PCPG (A)(i) Tile colour indicates the anatomical location of the tumour or, in the case of a metastasis, the associated primary tumour.An asterisk indicates whether the tumour had a non-chromaffin expression signature.(ii) Tile colour indicates the clinical behaviour of the tumour.Primary tumours resected from patients with metastatic disease are further annotated as to whether the tumour was able (filled circle) or unable (open circle) to be confirmed as the metastatic clone through sequencing of a paired metastasis.(iii) An identifier linking samples from the same patient.(iv) Indicates if the sample was resected after a cytotoxic treatment regime.(v) Immunohistochemical scoring of Ki67 expression, tumours with greater than 3% positive cells are indicated by a black square.(vi) Tumour mutation burden (mutations per megabase) and (vii) number of structural variants are indicated using a heatmap.(viii-x) Signature analysis was performed using the MutationalSignatures package for R. SBS signatures ranked (top to bottom) based on the mean proportional contribution of variants.Output was filtered to retain signatures that contributed 15% of mutations and at least 500, 50, or 10 mutations for SBS, ID, and DBS signatures, respectively, in at least one tumour.Colour indicates the proportion of mutations from each class attributable to the respective signature for a given sample.(xi) The presence of a TERT or ATRX mutation detected by WGS.(B) Aggregate and per-sample copy number aberrations detected by WGS in sympathetic (rows 2 and 3) and parasympathetic (rows 5 and 6) PCPG.Regions that were detected as significantly recurrently deleted (downward arrow) or amplified (upwards arrow) are indicated (rows 1 and 4).Clinical behaviour and the presence of whole genome doubling are annotated to the left of per-sample data (rows 3 and 5), while the percentage of genome altered and number of structural variants observed in each tumour are annotated on the right.Losses and gains as well as percent genome altered are expressed relative to sample ploidy rather than diploid.(C) Genes that were detected as recurrently altered by WGS and predicted to be drivers by Cancer Genome Interpreter across the A5 cohort (top) and the respective number of mutations observed in two previously published datasets (bottom).(D) Schematic showing the DNA and protein position of two somatic mutations observed in EPAS1.

Figure 4 .
Figure 4. TERT promoter mutations and structural alterations (A) sn-ATAC-seq covering the TERT gene and promoter (inset) region.Data is shown for tumour cells from two ATRX mutant tumours, three TERT mutant tumours, and one tumour with a TERT structural variant.Normal chromaffin, adrenocortical, fibroblast, and endothelial cells are shown for contrast.(B)(i) Circos plot describing somatic alterations detected by WGS in E123-M1.Outer ring indicates chromosome, second outer ring marks SNVs (blue: C>A, black: C>G, red: C>T, grey: T>A, green: T>C, pink: T>G) and Indels (yellow: insertions, red: deletions).Third outer ring indicates total copy number, (green: >2, red: <2), the fourth outer ring indicates the copy number of the minor allele (blue: >1, orange: <1).The center circle displays structural variants.(ii) Copy number segmentation data from the AFF4 and TERT gene regions in E123-M1 showing multiple short segmental copy-number changes due to chromothripsis.The break points contributing to an AFF4-TERT fusion are marked with red dashed lines.(iii) Schematic of AFF4-TERT fusion detected by WTS.(C) TERT expression (y-axis, log2 CPM) versus gene copy number at the TERT locus (x-axis) (D) TERT expression (y-axis, log2 CPM) versus the methylation status (x-axis) of probes in the TERT promoter and TERT hypermethylated oncological regions (rows 1 and 3).Density plot of probe methylation betavalues (rows 2 and 4).Point and line colour indicate TERT/ATRX mutation status.

Figure 5 :
Figure 5: Loss of ATRX leads to telomere dysregulation (A) Genomic (top) and protein (bottom) position of alterations detected in the ATRX gene by WGS.(B) Ratio of tumour to normal telomere read content (y-axis) in relation to ATRX and TERT mutation status (x-axis).Presence (red) or absence (black) of C-circles is indicated by the dot colour.Paired tumours from the same patient are joined by a dashed line.A two-sided Student's t-test was used to test for differences.Values were averaged across paired tumours prior to statistical testing.(C) Telomere variant repeats (TVRs) of the type NNNGGG (and the reverse complement) were detected in WGS data from tumour and matched normal using TelomereHunter.Count values were normalized against intratelomeric read count (green) or total read count (blue).The tumour/normal ratio of normalized counts (y-axis) are shown with respect to the presence or absence of detected c-circles (xaxis).Only TVR sequences that were significantly different using a Students t-test (p<0.05)applied to values normalized to total reads are shown.TVRs that were significant after false-discovery correction (Benjamini & Hochberg, p<0.1) are indicated with a double-asterisk.P-values before (p) and after (pBH) correction are shown.Data points are coloured to indicate a tumour/normal

Figure 6 :
Figure 6: TERT/ATRX-alterations and their association with metastatic progression (A) A schematic illustrating the case categories with presence or absence of TERT/ATRX mutations in primary and paired metastatic tumours.(B) A polyclonal primary harbouring ALT and pTERT mutations giving rise to metastases (i) The centre panel shows the variant allele frequency (VAF, y-axis) for all somatic mutations detected in a metastatic primary (E143-P1) and the two paired metastases (E143-M1/2) from patient E143 (x-axis).The colour of each dot denotes the genomic copy number status in the region of each variant.The horizontal lines indicate the purity of each tumour and can be read as a proportion from the x-axis.Grey lines connect mutations shared across the paired tumours while a pTERT mutation is highlighted with a red line.The top panel shows the tumour to normal telomere content ratio for each tumour and the colour of each bar indicates the presence or absence of C-circles.The panel on the right shows the genomic copy number status along each chromosome (y-axis) for each tumour (x-axis).(ii) A schematic illustration of the clonal evolution of metastatic disease in patient E143.The cell colour indicates the presence of the ALT phenotype (green) or pTERT mutation (red).(C) The tumour mutation burden (mutations per megabase, y-axis) observed in each tumour with respect to TERT/ATRX gene mutation status (x-axis).A one-tailed Student's t-test was used.The y-axis has been truncated to accommodate an extreme outlier which was excluded during statistical testing.(D)

Figure 7 :
Figure 7: Differential gene expression between TERT and ATRX-altered and non-metastatic tumours Differential gene expression analysis was performed contrasting non-metastatic primary tumours with either (i) ATRX mutant tumours, (ii) TERT mutant tumours, or (iii) all metastasis and confirmed metastatic primaries.A fourth contrast between TERT and ATRX mutant tumours was also performed.(A) An upset plot showing the intersection of genes that were significant (adjusted p-value < 0.05, log-fold-change > 1) in each contrast.Bar colour indicates gene-ontology association for protein coding genes or gene biotype for non-protein-coding genes.(B) Differential gene-expression between non-metastatic primary and metastatic tumours showing fold change (log2, y-axis) versus adjusted p-value (-log10, x-axis).Genes that were also significant in non-metastatic primary vs TERT-altered and non-metastatic primary versus ATRX-altered contrasts are coloured according to gene ontology annotation.(C) Heatmap (centre panel) showing genes that were differentially expressed in both the ATRX vs nonmetastatic primary and ATRX vs TERT contrasts.Annotation bars on the left indicate whether the gene was also found in a differentially methylated region (ATRX vs non-metastatic primary), and the correlation of the expression of MKI67 to the expression of each gene (Spearman correlation).The right panel shows the gene expression determined by snRNA-seq in cells aggregated by cell type (right sub-panel) or ATRX/TERT mutation status (left sub-panel).Dot colour indicates mean expression while dot size indicates the fraction of cells expressing the gene.(D) Expression (log2 CPM, y-axis) of genes found to be differentially expressed between ATRX-altered and non-metastatic primary tumours.Expression data is shown from the A5 cohort (top) and the TCGA/Flynn et al. cohorts (bottom).Data point colour indicates PCPG subtype.A Student's t-test was used to test for differences (E) Differentially expressed genes in both the TERT-altered vs non-metastatic primary and ATRX-altered vs TERTaltered contrasts.See panel C description for panel elements.The lower and upper hinges of each boxplot correspond to the first and third quartiles, respectively, and the median value is marked.The whiskers extend to the largest and smallest value no greater than 1.5 times the interquartile range above or below the upper and lower hinges, respectively.

Figure 8 :
Figure 8: Evolution of PCPG under treatment pressure (A) Clinical timeline of for patient E169.(Bi) Shared and private variants (y-axis) between paired metastases (x-axis) taken before and after CVD treatment.(Bii) Total number of variants detected by WGS (top left, bottom right) and number of variants shared between tumours (bottom left).(Biii) Shared and private structural variants between paired metastases.(Biv) Copy number status along each chromosome (y-axis) in paired metastases (x-axis).(C) Mutation signature analysis using COSMIC v3 SBS (i) and InDel (ii) signatures (y-axis) in paired metastases (x-axis).Heatmap colour indicates signature contribution.(D) Expression of MGMT in the A5 cohort.(E-G) Patient E167, see description for panels B-C.(H) Expression of MLH1 in the A5 cohort.(I) Total mutation counts (y-axis) for PCPG tumours (x-axis) in the Project GENIE data registry.Bar colour indicates the presence of a mutation in the mismatch repair pathway.(J) Trinucleotide context for mutations observed in E167-M1 (top), E167-M2 (second from top), the highest mutation load tumour from the GENIE dataset (second from bottom), and ctDNA derived from a patient with metastatic SHDB-related PGL treated with Temozolomide.(K) 18 F-FDG-PET imaging for Temozolomide treated patient at time of blood draw for cfDNA analysis.(L) Ichor CNA analysis of ctDNA derived from Temozolomide treated patient (M) Variant allele frequencies (y-axis) for somatic variants (x-axis) observed in ctDNA derived from Temozolomide treated patient.Datapoints are coloured to indicate transition/transversion or insertion/deletion type.